Steady-state spectra, current and stability diagram of a quantum dot: 
a non-equilibrium Variational Cluster Approach 



O 



u 

■i— > 

a 

fl 
O 

o 



in 
o 

<N 



Martin NussQ Christoph Heil, Martin Ganahl, Michael Knap, 
Hans Gerd Evertz, Enrico Arrigoni, and Wolfgang von der Linden 

Institute of Theoretical and Computational Physics, 
Graz University of Technology, 8010 Graz, Austria 
(Dated: July 25, 2012) 

We calculate steady-state properties of a strongly correlated quantum dot under voltage bias by 
means of non-equilibrium Cluster Perturbation Theory and the non-equilibrium Variational Cluster 
Approach, respectively. Results for the steady-state current are benchmarked against data from 
accurate Matrix Product State based time evolution. We show that for low to medium interaction 
strength, non-equilibrium Cluster Perturbation Theory already yields good results, while for higher 
interaction strength the self- consistent feedback of the non-equilibrium Variational Cluster Approach 
significantly enhances the accuracy. We report the current-voltage characteristics for different in- 
teraction strengths. Furthermore we investigate the non-equilibrium local density of states of the 
quantum dot and illustrate that within the variational approach a linear splitting and broadening of 
the Kondo resonance is predicted which depends on interaction strength. Calculations with applied 
gate voltage, away from particle hole symmetry, reveal that the maximum current is reached at the 
crossover from the Kondo regime to the doubly-occupied or empty quantum dot. Obtained stability 
diagrams compare very well to recent experimental data Phys. Rev. B, 84, 245316 (2011) . 
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I. INTRODUCTION 



The understanding of non-equilibrium phenomena 
in strongly correlated many-body systems may reveal 
previously unknown fundamental aspects of physics as 
well as prove crucial for the development of technical 
applications in the fields of nano or molecular electronics. 
Currently, combined insight of experiments on nano- 
devicesi — and results from artificial quantum simulators 
(see e.g. refs. 0-HI) are capable of providing coherent 
insight into the non-equilibrium behavior of the quantum 
world. These experiments provide both a challenge for 
theoretical concepts as well as an accurate check for 
theoretical results. Both nano-devices and condensed 
matter simulators are often described remarkably well 
by interacting model Hamiltonians which are in general 
not exactly solvable. 

Here we focus on a model of a single quantum dot, the 
single-impurity Anderson model (SIAM) 9 . This model, 
incorporating spin and charge fluctuations as well as 
Kondo correlations, has been studied as an idealized 
realization of an interacting system by a wide array 
of techniques in equilibrium (for an overview see e.g. 
Hewson ref . [To|) • 

However, the evaluation of dynamic quantities of 
strongly correlated quantum many-body systems out 
of equilibrium poses a notoriously difficult problem. A 
particular challenge for the SIAM is that it is expected 
to remain in a strong coupling regime, even under the 
influence of a bias voltageii. Additional challenges are 
posed by the need for particle and/or energy dissipation 
mechanisms. Techniques which allow the calculation of 
physical quantities beyond the linear response regime are 
quite restrictive up to now and no full understanding of 



the non-equilibrium dynamics or the steady-state under 
bias are available for the SIAM, apart from some special 
cases^r— . Nevertheless, controlled results have been 
obtained with the analytical Bethe Ansatz for some phys- 
ical quantities in the interacting resonant level modeU 5 - 
and for the SIAM±£. Perturbative calculations^ 7 -^ 
as well as non-crossing approximation studie s 19 ! 20 
extend early fundamental work on the non-equilibrium 
problem 2 ^ — and the impurity out of equilibrium^ 5 - — . 
Depending on the setup of parameters and model details, 
insight may be gained by semi-classical methods 2 ^ or 
master equation approaches^ 9 -. Recently moreover, 
several techniques, which have proven very successful 
in the equilibrium theory, have been extended to the 
non-equilibrium case. Among them are many-body 
cluster method s 30 ! 31 , renormalization group (RG) ap- 
proaches 3 -^—, flow equation method s 39 ! 40 , real time path 
integral calculations 4 ^, diagrammatic quantum Monte 
Carlo (QMC) 4 ^— or QMC methods based on a complex 
chemical potential 4 ^"—. The Gutzwiller approximation 
has been generalized to the time-dependent case 49 and 
so has numerical renormalization group (NRG) 5 -^— 
where however some issues with the use of Wilson 
chains in non-equilibrium systems have been pointed 
out by Rosch 54 . Dual fermion approaches 55 have been 
proposed as well as super operator technique s 56 ! 57 . 
Some recent work attempts to compare several of these 
theories 5 -^ — and shed light on the critical issue of time 
scales involved 6 ^. Finally, some results for the SIAM 
are available^ 2 - from numerically exact time evolution 
by a combination of Density Matrix Renormalization 
Group (DMRG) 63 and successive time evolution via 
time-dependent DMRG (tDMRG)^* 4 ^— . They are 
currently limited to small bias voltages and moderate 
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interaction strengths. 

In the present paper we explore the capabilities of 
non-equilibrium Cluster Perturbation Theory (nCPT) 30 
and benchmark the non-equilibrium Variational Cluster 
Approach (nVCA) 31 on the SIAM. We obtain the full 
current-voltage characteristics which we compare to 
results from a very accurate time evolution by means 
of Time Evolving Block Decimation (TEBD) 67 and 
QMC— . We also comment on the comparison to other 
results obtained by tDMRG^, perturbation theory-^ 
and FRG 33 . Because of good agreement with these 
data, we then proceed to evaluate the single-particle 
spectrum of the quantum dot in the steady-state. Since 
this is a dynamic quantity, it is even harder to obtain 
for most numerical methods. We find a linear splitting 
of the Kondo resonance which depends on the inter- 
action strength. Detailed results for the particle-hole 
symmetric model where Kondo correlations dominate 
are presented and supplemented by data in the whole 
parameter space including an applied gate voltage. We 
highlight the crucial edge which nVCA gains over nCPT 
via its self-consistent feedback. 

Both many-body cluster techniques, nCPT and nVCA 
are based on the well established equilibrium coun- 
terparts Cluster Perturbation Theory (CPTjSSiZtt and 
the Variational Cluster Approach (VCA)£!^. They 
make use of the non-equilibrium Keldysh-Schwinger 
Green's function technique^—. The present paper 
extends and benchmarks the ideas for adapting VCA 
to non-equilibrium situations introduced in ref. [3l|. We 
attempt to give a comprehensive overview of the current 
capabilities and shortcomings of nCPT and nVCA for 
the application to steady-state problems of strongly 
correlated systems. The SIAM provides an excellent 
probing ground for our purposes as a model where the 
effects of correlations are crucial. It is however still 
relatively simple which permits systematic study and 
some results are available which allow for comparison. 
Reasonable results for the SIAM in equilibrium have 
been obtained previously by CPT as well as VCA—. 
Both cluster methods are approximate but they yield 
fairly reliable results and are therefore interesting due 
to their great flexibility and versatility. They are com- 
putationally not very demanding and allow in principle 
to treat a wide range of fermionic / bosonic lattice 
Hamiltonians out of equilibrium. Possible extensions 
include electronic multi-band or multi-orbital systems in 
one, two or three dimensions also including phonons. 
This paper is organized as follows: In sec. |TT] we sketch 
the SIAM and describe the setup used. We proceed 
by outlining nCPT (sec. UTTAj) and nVCA (sec. ITTTB]) . 
In sec. IIV A[ a comparison of steady-state currents 
as obtained by nCPT, nVCA, DMRG and successive 
TEBD 77 and QMC^ is presented. The non-equilibrium 
local density of states (nLDOS) is examined in sec. UVBl 
Finally effects of an applied gate voltage on the steady- 
state current are studied in sec. IIVCI 



II. MODEL OF A SINGLE QUANTUM DOT 

We model the setup, consisting of a single correlated 
quantum dot in-between two metallic leads, by a single- 
site Hubbard model embedded in a one dimensional infi- 
nite tight-binding chain. The Hamiltonian of this single- 
impurity Anderson model (SIAM}2 reads (see fig. [J) 

T~L — ^dot + Hlead + H CO up 

^dot = e/ ftf* + Un{h{ (lb) 



^lead — I €a Z-^ ° iacr Cia(J ^ Ciacr C J a(T I 

a,a \ i=0 (ij) J 

(lc) 

Ticoup = -f £ (c\ aa f a + fl <w) • (Id) 

a, a 

The electronic annihilation (creation) operators c iaa ,f a 
(c\ aa , fl) obey the usual anti-commutation relations with 
spin a = {t, 1} and annihilate (create) electrons in the 
left or right lead a = {L, R} or the quantum dot, re- 
spectively. The particle number operator of the quan- 
tum dot is defined by = f£f a . U represents the 
on-site Hubbard repulsion. The single-particle energy 
of the quantum dot €f = —^ + Vq is comprised of a 
gate voltage Vq and a shift by — ^ (eq. (|lb[) ). such that 
the particle-hole symmetric point is reached when the 
gate voltage vanishes and the bias is applied in an anti- 
symmetric fashion as described below. Non-interacting 
left and right leads are described by tight binding chains 
with a nearest-neighbor hopping t (eq. (fTc|) ). The lead 
on-site energies are denoted by e^/^. A bias voltage Vb 
may be applied to the system in an anti-symmetric fash- 
ion by setting cl = — —£r — ~Hr — -^F, where 
Hl/r denote the chemical potentials of the decoupled 
leads (t f = 0). Finally, the symmetric tunneling be- 
tween the non-interacting leads and the quantum dot is 
denoted by t' (eq. (|ld[) ). In this work we adopt units in 
which h, e, as well as the inter- lead hopping t are equal 
to 1. The lead-dot tunneling is fixed to t' = 0.3162 
throughout this paper, which implies an Anderson width 
of A = 7rt /2 piead(0) = ^ « 0.1. The prefactor for the 

current is then given by jo = In the following, we 

always consider the zero temperature case. 



III. METHODS 

Upon introducing the basic notation we follow stan- 
dard literature (see e.g. ref. [T8|) . The non-equilibrium 
single-particle properties are provided by the single- 
particle Green's function in the 2x2 Keldysh space (de- 
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FIG. 1: (Color online) Illustration of the SIAM as prepared 
for use within nCPT and nVCA. An interacting quantum dot 
(QD) is located in the middle of two non-interacting electronic 
leads. In order to employ the cluster approaches, the system 
is divided into three pieces, a left part, a right part and a 
central region of length L, which includes the quantum dot 
as well as parts of left and right leads. 



noted by tildes) 



G 



G* G K 
G A 



(2) 



where 

qR/A/K 

are again matrices in site/spin space and 
functions of two time coordinates^. R denotes the re- 
tarded, A the advanced and K the Keldysh component 
of the single-particle Green's function which for (non- 
relativistic) fermions read 

G^(r 1 -r 2 ) = z^(r 1 -r 2 )([q(r 1 ), C t(r 2 )]) 



r 2 ) = -ze(r 2 -r 1 )([c l (r 1 ),c](r 2 )}) 
r 2 ) = -licltn^fa) + c^T^fa)) , 



where z, j denote site as well as spin, ti,t 2 real time and 
[A, B] = AB — BA the standard commutator. Note that 
in the following we consider the spin symmetric model 
(see eq. ([Taj)) and therefore suppress spin-indices. In the 
formulas for the current we therefore include an addi- 
tional factor of two. The handy matrix relations G a (ti — 
T2 ) = (G fl )tfo-Ti) and G K (n-r 2 ) = -(G if )t(r 2 - n ) 
follow from the definitions. 

Our goal is to investigate the steady-state in which the 
system becomes time-translationally invariant. There- 
fore we may Fourier-transform to the energy domain and 
evaluate single-particle steady-state expectation values 

The steady-state current through site m can be evalu- 
ated^^ by the time derivative of the total particle num- 
ber to the left of this site 



m—l 



j = (N M (t)) , N M 



= E 



This expression for the current may be rewritten in terms 
of the Keldysh component of the single-particle Green's 
function 



Ji i+1 — ti i+1 



duo 
2tt 



5Re(G^ +1 H-Gf +H H) , (3) 



where ta+i is assumed to be real. The non-equilibrium 
density of states at site i is given by 

7T 

which has to fulfill the spectral sum rule 



1 = 



Pi(uj) Mi. 



(4) 



Since all this information is contained in a single object, 
the single-particle Green's function in Keldysh space G 
(eq. (j2j)), we seek a way of evaluating it. This is in gen- 
eral impossible to do exactly but a handle on G is given 
by nCPT and nVCA. These approximations are to be 
discussed in the following. 



A. Non-equilibrium cluster perturbation theory 

The idea of nCPT is to split an (infinitely) large system 
T~L for times r < To into a set of small decoupled clusters 
(described by the Hamiltonian h), for which the single- 
particle Green's functions can be determined exactly, by 
numerical means 

H = h + 6(T-T )f. 

At time To the coupling of the individual subsystems 
T = % — h is switched on and the solution for the single- 
particle Green's function of the total system can be ob- 
tained by the CPT equation 69 in Keldysh-, site- and spin- 
space 



G 



n. 



(5) 



where G/g denote the single-particle Green's function in 
Keldysh space of the total system/split systems. The 
unit matrix in Keldysh space is denoted by IL. Eq. (|5j) 
is a strong coupling perturbation theory result and holds 
up to first order in the inter-cluster hopping T, as far 
as the self-energy is concerned. As a consequence, the 
self-energy of the total system is approximated by 
the self-energy of the initial system E^. This implies 
that nCPT/nVCA become exact in the non- interacting 
limit, since in this case, the self-energy of the total system 
agrees with that of the initial system, which are both 
zero. The results of nCPT and nVCA converge towards 
the exact results with increasing cluster size. 
In the present paper we illustrate this approach for the 
SIAM out of equilibrium. We start out by splitting the 



infinite chain (eq. ([Taj)) into three parts: i) an interacting 
central region of length L consisting of the interacting site 
(the quantum dot) as well as an equal amount of sites 
of the left and the right leads, ii) a semi-infinite, non- 
interacting left region consisting of the remaining part 
of the left lead and iii) a semi-infinite, non-interacting 
right region consisting of the remaining part of the right 
lead (see fig.[Tj) . To proceed, the single-particle Green's 
functions in Keldysh space have to be determined exactly 
for those three systems. Results for the left and right 
part are available analytically by the retarded Green's 
function of a semi-infinite tight binding chain^l 

sS^m = v ii R j(u) - v o!i+j(u) ( 6 ) 

L/R, v -I ( UJ-6 L / R 

^ M = vw-(.-ev* )2 { — wr (T) 

where Vij is the retarded Green's function of the infinite 
tight binding chain. The single-particle Green's function 
of the central interacting region can be determined in 
the Q-matrix formalism by the Band-Lanczos algorithm 
(see ref. [76| and ref. l82l ). The advanced component can 
always be obtained by the relation G A = (G^)^. Before 
coupling the three subsystems at time r = To, each of 
them is in equilibrium at different chemical potentials 
Vl/r/c- Therefore we may evaluate the corresponding 
Keldysh components by the relation^ 

G*(w, (j) = (G r (cj) - G a (uj)) (1 - 2p FD ( w , (j,)) , (8) 

where pfd denotes the Fermi-Dirac distribution function 
at inverse temperature f3: pfd = 1Jre ^-^ • In the zero 
temperature limit, (1 — 2£>fd(<^, aO) may be re-expressed 
as sign(o; — /x). This is the only expression in which 
the chemical potential /i enters. It is crucial that this 
relation does not hold in a non-equilibrium situation 
any longer. The hermitian part of G K (o;,/i) is zero in 
equilibrium. The imaginary part consists of contribu- 
tions due to delta peaks for finite size systems like the 
central regions. The operator T just contains the two 
hopping terms from the left to the central and from the 
central to the right region. At time To, the hopping 
processes between these three regions are switched on 
and the steady-state single-particle Green's function G 
is determined using nCPT, eq. ([5]). As in the equilibrium 
case, the CPT results can be improved by the variational 
cluster approach, in which the single-particle part of the 
initial system is suitably modified. In the following a 
variational extension (nVCA) of the scheme described 
above will be presented following ref. |3~0 . 



B. Non-equilibrium variational cluster approach 

The idea of enlarging the size of the central region in 
nCPT is to find a better approximation for the starting 
point of perturbation theory. As pointed out before, in 
CPT the self-energy of the total system is approximated 
by that of the decoupled system. By increasing the size 
of the central region, the self-energy of the decoupled sys- 
tem converges gradually towards that of the total system. 
In nVCA an even more suitable initial system is chosen. 
This can be achieved by making use of the fact that the 
decomposition of T-L into clusters h and inter-cluster parts 
T is not unique. We are free to add single-particle opera- 
tors A(x), which depend linearly on parameters x, to the 
initial hamiltonian h, provided we subtract them again 
from f 

H=(h + A(x)) + 0(t - r ) (f - A(x)) . 

Here, 

A(x) = x i , 
i 

where is quadratic in the fermion operators. This 
parametrized single-particle field introduces additional 
freedom in the starting guess for perturbation theory and 
can be used for a self-consistent feedback on the clusters. 
The system described by (ft + A(x)^ , usually referred to 

as reference system in the context of VCA, has the same 
structure as in equilibrium VCA. The condition to fix 
the single-particle parameters x, however, is different in 
nVCA and VCA. Equilibrium VCA is based on the Self- 
energy Functional Approach (SFA) 84 i 85 and provides a 
variational principle for the generalized grand potential 
functional, which is not well defined in the context of non- 
equilibrium systems any longer. An alternative criterion 
for fixing the variational parameters x was introduced in 
ref. [3l| and compared to the traditional VCA criterion in 
ref. [76|. It is based on the idea of starting from a system 
which is as similar as possible to the total, original system 
in terms of physically observable quantities. We demand 
the expectation values of the operators Ai to coincide in 
the initial (reference) system and the steady-state of the 
total system, i.e. 

<A,) g =L (A,) G . 

For example, adding variational freedom in the on-site 
energy of the quantum dot, corresponding to A = x e/ n£ 

yields the self-consistency condition: (n£) g = 
From a more conceptual point of view, it is interesting 
that these self-consistency conditions follow from the con- 
dition^ 
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which is closely related to the stationarity condition of 
the generalized grand potential functional Q in the equi- 
librium approach 86 . Here t\ = (^qJ * s a P aim matrix 

in Keldysh space and the subscript zero denotes non- 
interacting Green's functions. From the numerical point 
of view, one has to find the roots of an n-dimensional set 
of non-linear equations (where n is the number of varia- 
tional parameters x). 

Like in first order Dirac perturbation theory, the influ- 
ence of the perturbation increases with time (perturba- 
tions introduced due to the coupling with hopping ele- 
ments t at To are proportional to t • (r — tq) since we 
are considering first order perturbation theory) and one 
might argue that nCPT is bound to fail in the long-time, 
steady-state limit, even for small couplings between the 
clusters^. This argument cannot be true in general, as 
nCPT yields exact results in the case of non-interacting 
particles, although the initial decoupled systems are far 
from the steady-state behavior. 

The initial system in nCPT is independent of the non- 
equilibrium situation and this shortcoming is improved to 
some extent within nVCA, where the information about 
the applied bias voltage is self-consistently fed back to 
the initial reference system. As we will see in this paper, 
there are circumstances under which nCPT already yields 
reasonable results. In general, however, we observe that 
nVCA represents a significant improvement over nCPT. 
This implies, on the one hand, that under steady-state 
conditions, the self-energy in the central cluster is signif- 
icantly modified when going from nCPT to nVCA. On 
the other hand, it indicates that the steady-state situ- 
ation can be mimicked in an equilibrium system by the 
auxiliary one-particle terms, determined self-consistently. 
A point for improvement of this approach is to modify 
the self-energy so that it is a genuine non-equilibrium 
one. 

In this paper we compare nCPT with nVCA for different 
auxiliary one particle terms A(x). In agreement with our 
previous experience, that the variational hopping param- 
eters are crucial for the physics in the Kondo regime^ 6 -, 
we observe that it suffices to consider only hopping pro- 
cesses inside the central region of the reference system. 
The most general form of the auxiliary one-particle term 
considered in this work reads 

A(x t ,x t /,x t J = x t > ( c Lr fa + he) 

a, a 

+ X * Yl C i-laa + hc ) 

+*t b J2(4 aa c e _ laa +hc) , 

where £ = (L — 3)/2. It contains three different hopping 
processes. The first term describes the hopping processes 
to and out of the quantum dot. The second term contains 
nearest neighbor hopping processes between the lead sites 



of the central region except the border one. The last hop- 
ping process involves the nearest neighbor sites at the 
boundary of the central region. In the sequel, we con- 
sider the following cases: (uVCAt) with x t ' = x t = x tb , 
(nVCA tjt /) with x tb = 0, (nVCA t /) with x t = x tb = 0, 
and (nVCAt b ) with Xt = Xf = 0. 

A bias voltage Vb is applied in an anti-symmetric man- 
ner by setting the on-site energies as well as the chem- 
ical potentials of the left and right leads to cl = [II — 
—6r = —fiR = -^p-. Note that the sites of the leads 
which are incorporated in the central region of the refer- 
ence system are also subjected to these on-site energies. 
All calculations are done for sizes of the central region 
of L = 3, 7 and 11-sites, which corresponds to symmet- 
ric central regions. We note that L = 5, 9 suffer from 
a finite-size gap which closes with increasing L. This 
gap arises due to an even amount of sites to the left as 
well as to the right of the quantum dot. We do not con- 
sider even L since this would correspond to an unequal 
amount of sites of the left and right lead in the central 
region. Thereby, the geometric symmetry of the problem 
would be spoiled and the application of a bias voltage 
which respects particle-hole symmetry would not be pos- 
sible. 

In order to reduce finite size effects of the calculated 
quantities and to make an extrapolation to infinite cen- 
tral cluster sizes easier we scrutinized the idea of aver- 
aging over boundary conditions as outlined in ref. [H3 or 
ref. [88|. In CPT/VCA one has the freedom to add single- 
particle terms to h and subtract them again in T. There- 
fore we may add hopping terms which connect the first 
and last site of the central cluster and which have an 
arbitrary phase. The spectral function of the impurity 
site for a central region with complex periodic boundary 
conditions is in general no longer particle hole symmet- 
ric, which makes it necessary to average G^(cj) over <p to 
retain the symmetry. 

The result, however, was disappointing as the expected 
faster convergence towards the thermodynamic limit has 
not been confirmed. 

In this work we do not consider the long-range part of the 
Coulomb interaction which leads to additional charging 
effects in real devices. 



IV. RESULTS 

In the following we compare nCPT and nVCA data 
for the steady-state current with TEBD^i and QMG^ 
results in the particle-hole symmetric model. We 
elaborate on the non-equilibrium density of states and 
finally discuss results for the steady-state current away 
from particle- hole symmetry by applying a gate voltage. 
Earlier preliminary results obtained for the SIAM out of 
equilibrium by nCPT and nVCA are available in ref. [89| 
and ref. [9(1 
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FIG. 2: (Color online) Steady-state current- volt age characteristics in the particle-hole symmetric case (Vb = 0). Results are 
shown for small to large interaction strength: U = 4 A (top row), U = 8 A (second row), U = 12 A (third row) and U = 20 A 
(bottom row). The first column contains the full current- volt age characteristics as obtained with nCPT for L — 3, 7 and 11. As 
benchmarks, accurate TEBD results^ and the linear response result jii n — 2GqVb are depicted as well. In some bias regions, 
only an upper bound can be obtained by TEBD. This upper limit is depicted as a magenta bar. In the second column a close-up 
of the low bias region is presented. Here we compare nCPT (L = 3, 7 and 11), nVCAr (L = 3, 7 and 11) with TEBD data as 
well as QMC results^. Note that the QMC data have been obtained in a wide-band limit. The third column contains results 
obtained for L = 7 for various sets of nVCA variational parameters: nCPT, nVCAT, nVCA t /, nVCA t t /, nVCA^. Legends are 
displayed once in the top row for the respective column. 
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A. Steady-state current 

Here we investigate the steady-state current-voltage 
characteristics of the particle-hole symmetric SIAM. In 
this parameter region, Kondo correlations are especially 
important. All nCPT and nVCA results are for an in- 
finite system using a self-energy based on an L < 11 
site interacting reference system. Although the length 
of Kondo correlations scales exponentially in interaction 
strength, it has been shown in ref. [76| that the VCA 
approximation, similarly to other approximate methods 
such as FRG or Gutz wilier wave functions, is capable 
of retaining qualitatively most features of the ongoing 
Kondo physics, although, possibly, with renormalized 
scales. The steady-state current can be evaluated on any 
link either within the central region or between the cen- 
tral and the neighboring regions, since we find that the 
continuity equation is fulfilled within at least 10 -6 rela- 
tive to the steady-state current amplitude. We note again 
that in the non-interacting case nCPT as well as nVCA 
become exact so we do not show these data explicitly. 

Results obtained by eq. (|3j) for the nCPT case are 
shown in the left column of fig. [2] for various values of 
the interaction strength U. The linear response result 
jiin = 2 Go Vb , where Go is the conductance quantum, is 
fulfilled within nCPT (and nVCA) in the very low bias re- 
gion. This is because these methods fulfill the Friedel sum 
rule in equilibrium^. Actually, nCPT tends to stay in 
the linear response regime longer than nVCA and there- 
fore overestimates the current for low (but not very low) 
bias voltages. We model the leads by one dimensional 
tight binding chains, yielding a semi-circular density of 
states (eq. Q) of bandwidth D = 40 A. This implies a 
vanishing current in the high bias limit at Vb = 40 A, 
due to non-overlapping lead density of states. This limit 
again is trivially fulfilled within nCPT and nVCA. 
For comparison, numerically exact results 67 as obtained 
by a real-time evolution with TEBD are also indicated in 
fig. El Note that the TEBD data are obtained from the 
steady-state plateau in the time dependence of the cur- 
rent. At small and medium bias, the current converges 
well and the TEBD results provide an accurate bench- 
mark. At higher bias, convergence is low , but the TEBD 
data still provide a reasonably reliable upper bound for 
the steady-state current. 

In the intermediate bias region it is interesting to in- 
vestigate the behavior of nCPT with increasing size of 
the central region L. As can be seen from the plots, 
for any interaction strength increasing L yields mono- 
tonically improving results. While for low interaction 
strength U = 4 A the nCPT results almost coincide with 
the TEBD data, greater deviations arise at higher inter- 
action strengths. For very large interaction strength (see 
for example U = 20 A), the lengths of the central region 
L considered here are not sufficient. At high bias voltage 
some spurious finite size effects of the reference system 
are visible in the form of peaks in the steady-state cur- 
rent. 



Next we would like to discuss the performance of 
nVCA. An illustration is shown in the mid column 
of fig. EJ Here we compare nCPT and nVCAT for 
L = 3, 7 and 11 again for increasing interaction strengths 
U /A = 4, 8, 12 and 20. We show a zoom to the low bias 
region where benchmarking data from other techniques, 
like TEBD data 6 - 7 -, tDMRG^, FRG^, and QMC data 
obtained by Werner et al 42 in the wide-band limit are 
available. One aspect to note immediately is that nVCA, 
like TEBD, departs sooner than nCPT from the linear 
response data, which is due to the better reproduction 
of the Kondo resonance within nVCA. The thinning of 
this narrow resonance at zero bias with increasing in- 
teraction strength is responsible for the departure from 
the linear response result—. Since the Kondo resonance 
is better accounted for in nVCA, the curves leave the 
linear response data sooner than the respective nCPT 
results. The TEBD data 67 and QMC data obtained by 
Werner et al^ in the wide-band limit are also plotted 
and serve as a benchmark up to Vb ~ 5 A before effects 
of a different lead density of states become important. 
Curves obtained by tDMRG^ / FRG^ are not shown 
but lie basically on top of the TEBD / QMC data. Early 
data from perturbation theory^ 7 - are known to have some 
additional spurious bumps in the low bias region. The 
improvement due to the variational feedback in nVCA 
becomes clear by comparing to the nCPT results for the 
same L, especially at higher interaction strengths. 
The right column of fig.[2]shows a comparison of the per- 
formance of different variational parameters considered 
within nVCA for L = 7. For small interaction strength 
all nVCA parameter sets work equally well. With in- 
creasing interaction strength, however, the different vari- 
ational parameters predict a different behavior of the 
steady-state current. It should be noted here that some 
parameters like the ones used in nVCA tb cause almost 
no deviation from the nCPT result while others improve 
it appreciably, like e.g. nVCAT- It predicts a two peak 
structure in the high bias voltage regime, which how- 
ever is not observed directly in the TEBD benchmark 
data. No final conclusion can be drawn about the behav- 
ior of the current in this bias regime because TEBD can 
only predict upper bounds for the steady-state current 
there. The position of the dip in the steady-state current 
is at Vb = U for all interaction strengths. We note that 
this is where the bare level position of the quantum dot 
€f = — is about to stop overlapping with one of the 
lead density of states (while still overlapping with the 
other). 

The calculation for two independent variational param- 
eters (nVCAt 5 £/) yields similar results as nVCAx in the 
lower bias region but goes to zero quickly for high bias 
voltages and does not show a dip. 

On the whole, one may say that the self-consistent feed- 
back implemented within nVCA provides a significant 
improvement over the nCPT results. This has been al- 
ready observed in ref. 76. However, there this result was 
deduced from the convergence of results with increas- 
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FIG. 3: (Color online) Non-equilibrium local density of states of the quantum dot pf as a function of energy uj and bias voltage 
Vb in the particle- hole symmetric case (Vg = 0). Region A represents the experimentally most interesting region where no 
effects of the lead bands come into play and one observes a splitting of the equilibrium Kondo resonance. Region B is dominated 
by high energy incoherent excitations which are expected to form continuous Hubbard bands in the thermodynamic limit. In 
region C, the lead density of states becomes relevant, causing a new low energy resonance for L >= 7. Data are depicted for 
U = 12 A and an artificial numerical broadening of + = 0.05. nCPT results are shown in the first row for system sizes of L — 3 
(left), L = 7 (mid) and L = 11 (right). The second row displays the corresponding nVCAr data. For reasons of visualization 
the color-scale is cutoff at p ma x = 2. Note that a horizontal line through the center of each plot at V = corresponds to the 
equilibrium density of states. The spectral function sum rule eq. © is fulfilled for each applied bias voltage. The L = 11 results 
are noisy as fewer data points have been obtained. 




co/A 



FIG. 4: (Color online) Non-equilibrium local density of states 
of the interacting quantum dot. Data are shown for U — 12 A 
and an artificial numerical broadening of + = 0.05. Results 
are obtained by nVCAr with L = 3. 



ing cluster size and not with a benchmark comparison 
with alternative numerical methods. For low interac- 
tion strengths U < 4 A, bare nCPT already performs 
very well (independent of L) while for larger interaction 
strengths the variational improvement of nVCA becomes 
important. Motivated by the success of nVCA for the 
steady-state current, we proceed by evaluating the non- 
equilibrium local density of states of the quantum dot. 



B. Non-equilibrium local density of states 

The nLDOS in the quantum dot pf(oo) as obtained by 
nCPT and nVCAr is depicted in fig. [3] for three sizes of 
the central part of the reference system: L = 3, 7 and 11 
and a large interaction strength of U/A = 12 (corre- 
sponding to the steady-state current in fig. [2] (third row)). 
We plot the nLDOS in a density plot as a function of en- 
ergy uj (horizontal) and applied bias voltage (vertical). 
The equilibrium result, consisting of a thin Kondo reso- 
nance at the Fermi energy (ep = here) and two broad 
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FIG. 5: (Color online) Interaction dependent splitting of the 
Kondo resonance at uor under bias. Data shown are for differ- 
ent interaction strengths U = 4, 8, 12 and 20 A (solid lines). 
The resonance merges with the incoherent high energy spec- 
trum (Hubbard bands), located at ujh ~ =b-j (dotted lines) 
at a certain bias voltage V m ~ 15, 17, 21 and 28 respectively. 
These values have been obtained by the linear fits (indicated 
as dashed lines). If the Kondo resonance would pin at the 
chemical potential of the leads cjk, would follow the linear 
behavior indicated as /jll/r (dash dotted line). 



incoherent peaks located at « =b^ with a width of ~ 2 A, 
can be inferred from a horizontal cut at Vb =0. For fi- 
nite bias voltage, a splitting of the Kondo resonance is 
observed in both nCPT (top row) and nVCAr (bottom 
row). It is well known that the noncrossing approxima- 
tion (NCA) predicts a splitting of the Kondo resonance 
into two under voltage bias and that within second order 
perturbation theory the resonance is not split but sup- 
pressed only*i£ A linear splitting and slight broadening 
of the Kondo resonance with increasing bias voltage is 
proposed e.g. in ref. [9l|. Intuitively it is expected that 
the split Kondo resonances pins at the chemical poten- 
tials of the leads. Several other methods yield a split- 
ting with different features: real time diagrammatic 92 
and scaling methods 93 as well as the equation of motion 
techniqu o 19 i 2Q i 94 ~ — . Within fourth order perturbation 
theory the Kondo resonance splits into two, which are 
located near the chemical potentials of the two leads^. 
nVCA predicts a linear - interaction dependent splitting 
(see fig. H]), while nCPT yields a splitting which shows 
a roughly quadratic dependence on Vb- In contrast to 
the prediction of ref. 91 we do however not observe a 
simple pinning of the Kondo resonance at the chemical 
potentials of the leads. Our results indicate that the po- 
sition ujk of the split Kondo resonance depends on the 
interaction strength U: ujk = ±S^Vb for Vb < V m 
and ujk = ±Tf for Vb > V m (see fig.EJ). Here V m is 
the voltage where the Kondo resonance merges with the 
high energy part of the spectrum located at ujh ~ ±¥: 



Vm = 15, 17, 21 and 28 A for U = 4, 8, 12 and 20 A. The 
U dependent values of Vm have been determined from 
the respective linear fit to the data. 
For high voltages, these split peaks merge with the Hub- 
bard bands and saturate, which has also been observed 
in fourth order perturbation theory calculations^. In 
this bias region a new low energy excitation is observed 
for L > 3 within nVCAx- This additional peak in the 
nLDOS has a dominant contribution to the total weight 
and is responsible for the two-peak structure observed in 
the nVCAT steady-state current (see fig.|5J). 



C. Finite gate voltage: steady-state current and 
stability diagram 

It is interesting to investigate the steady-state current 
away from the particle-hole symmetric point and the 
region where Kondo correlations are present. In this 
section the current through the quantum dot under bias 
is analyzed as a function of an applied gate voltage 
Vg and applied bias voltage Vb at fixed interaction 
strength U. Results are presented for L = 3. In this case 
no two-peak structure has been found in the nVCAT 
current in fig.EJ which is corroborated by fig. [3] (bottom 
left). The dependence of the current on the gate- 
and bias voltage as obtained by nCPT and nVCAT 
is depicted in fig. [6] for interaction strengths U/A = 4 
and U/A = 20. First we discuss the so called stability 
diagram (differential conductance as a function of bias- 
and gate-voltage) which is shown in the third column 
of fig. El The dashed squares mark the region which 
is typically accessed experimentally. In the nVCA^ 
result, the Kondo region (vertical line in the center) is 
reproduced extremely well, as are the Coulomb blockade 
regions (oc Vb) above and below. Our data compare 
very well to recent experimental data of ref. Q (figure 5 
therein). Effects at higher bias voltage arise from the 
finite bandwidth used here and are typically not seen in 
experiment. 

The steady-state current as a function of bias and 
gate voltage is shown in the first and second column 
of fig. [6] It is interesting to observe that the largest 
current is obtained exactly at the crossover points from 
the Kondo to the empty or doubly occupied quantum 
dot (these regions are marked by black-dashed lines in 
middle and right panel of fig. [6]). This aspect can be 
seen more clearly in fig. [3 where the current is plotted 
as a function of gate voltage. For U = 4 A there is 
not much difference between the nCPT and nVCA 
results as can be inferred from fig. [6] For U = 20 A, 
however, we see that the feedback mechanism in VCA 
has a significant impact and leads to smoother j — Vb 
curves due to suppression of finite size effects originating 
from the reference system. The sharp jumps at the 
crossover point in fig. [7] originate from abrupt changes of 
the particle number in the ground state of the central 
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FIG. 6: (Color online) Stability diagram and steady-state current as a function of bias voltage Vb and gate voltage Vg- In 
the first column, the steady-state current is shown as a function of bias voltage for various gate voltages. The second column 
contains a density plot of the steady-state current in the full Vb /Vg parameter space. In the third column a similar density 
plot is shown for the differential conductance G — (stability diagram). The four rows depict results obtained by nCPT 
for U = 4 A, nVCA T for U = 4 A, nCPT for U = 20 A, and nVCA T for U = 20 A. The squares in the two stability diagrams 
bottom right mark the experimentally most interesting region. In the nVCAT result, the Kondo region (vertical line in the 
center) is reproduced extremely well, as are the Coulomb blockade regions (oc Vb) above and below. For comparison to recent 
experimental data see ref. 2 (figure 5 therein). Effects at higher bias voltage arise from the finite bandwidth used here and are 
typically not be seen in experiment. 



11 



1.5 



CD 



=3 
O 



I 

t 






v R 

B 


=4A 




B 


= 16A 


1 


...v B 


=24A 


1 


V B 


=36A 








1 

i\ 
i 


M 
M 






\ ^ 






10 



20 



gate voltage V Q /A 



FIG. 7: (Color online) Dependence of the steady-state cur- 
rent on gate voltage Vg- The largest steady-state current is 
obtained in a parameter regime which crosses from the Kondo 
plateau to the doubly occupied or empty quantum dot. Data 
as obtained by nVCAT (L = 3) are presented for an interac- 
tion strength of U/A = 12. 



cluster L = 3. We expect that this step smoothens with 
increasing L. 

Concerning the reliability of the methods used, our data 
suggest that outside the parameter region where Kondo 
correlations dominate the single occupied quantum 
dot (in between the horizontal dashed black lines in 
fig. [6]) nCPT and nVCA perform almost equally well. 
Results are significantly easier to obtain outside the 
Kondo plateau (which was mentioned before e.g. in 
ref. [76| within exact diagonalization/CPT/VCA or in 
ref. [62| within DMRG). Thereby the convergence within 
system size is greatly enhanced with respect to the 
Kondo region and very accurate results may be obtained 
already with small systems. Therefore, we may argue 
that for the steady-state of the SIAM, nCPT and nVCA 
perform quite well outside the Kondo region for any 
interaction strength as well as in the Kondo region 
for small interaction strengths. On the other hand, in 
the Kondo region nVCA outperforms nCPT for higher 
interaction strengths. 



the thermodynamic limit which is necessary to account 
for particle and energy dissipation mechanisms. 
Results for the particle-hole symmetric model, which is 
dominated by Kondo correlations, have been compared 
to Time Evolving Block Decimation and quantum Monte 
Carlo. At low values of interaction strength they show 
excellent agreement already for non-equilibrium Cluster 
Perturbation Theory. For higher values of interaction 
strength, the self-consistency implemented within the 
non-equilibrium Variational Cluster Approach proves 
crucial in order to obtain reasonable results. Both 
methods coincide with the low bias linear response data 
for the steady-state current. Both methods furthermore 
become exact in the non-interacting limit. 
The non-equilibrium local density of states of the 
quantum dot exhibits a linear and interaction dependent 
splitting of the bias voltage within the non-equilibrium 
Variational Cluster Approach which is not visible in 
the non-equilibrium Cluster Perturbation Theory. At a 
certain (interaction dependent) bias voltage we find that 
this split Kondo resonance merges with the high energy 
incoherent part of the spectrum. 

When applying a gate voltage and thereby leaving the 
Kondo regime, calculations become a lot easier and 
non-equilibrium Cluster Perturbation Theory and the 
non-equilibrium Variational Cluster Approach appear 
to perform very well, which can be inferred from the 
convergence of our data. The highest current amplitude 
is obtained at the crossover from the Kondo to the 
empty or doubly-occupied quantum dot. Experimental 
stability diagrams are reproduced very well within the 
variational approach. They show a clear Kondo regime 
and a Coulomb blockade region. 

We may conclude that the non-equilibrium Varia- 
tional Cluster Approach is a promising method for 
the evaluation of steady-state quantities of strongly 
correlated model systems. It is both fast and versatile. 
Dynamic quantities become available in the whole 
complex plane and, in principle any fermionic or bosonic 
lattice model may be treated including multi-band or 
multi-orbital systems. Subjects for further investigation 
are the self-consistency criterion, which is not unique, 
as well as the use of different variational parameters 
and an optimized out of equilibrium variational principle. 



V. CONCLUSIONS 

We have presented results for the steady-state of 
the single-impurity Anderson model. We have applied 
non-equilibrium Cluster Perturbation Theory and its 
variational extension, the non-equilibrium Variational 
Cluster Approach to this model for a single quantum 
dot under bias. Both methods make use of the Keldysh 
Green's function formalism and are capable of working in 
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